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Abstract 

Green-Kubo and Einstein expressions for the transport coefficients of a fluid in a 
nonequilibrium steady state can be derived using the Fluctuation Theorem and by assuming 
the probability distribution of the time-averaged dissipative flux is Gaussian. These 
expressions are consistent with those obtained using linear response theory and are valid in 
the linear regime. It is shown that these expressions are however, not valid in the nonlinear 
regime where the fluid is driven far from equilibrium. We advance an argument for why 
these expressions are only valid in the linear response, zero field limit. 



I. INTRODUCTION 

In 1993 Evans, Cohen and Morriss [1], ECM2, gave a quite general formula for the 
logarithm of the probability ratio that in a nonequilibrium steady state, the time averaged 
dissipative flux takes on a value, J_,_(t), to minus that value, namely, J_(t) = -J+(t) . That is 
they gave a formula for ln[p(J_,_(t)) / p(J_(t))] from a natural invariant measure [1, 2]. This 
formula gives an analytic expression for the probability that, for a finite system and for a 
finite time, the dissipative flux flows in the reverse direction to that required by the Second 
Law of Thermodynamics. The formula has come to be known as the Fluctuation Theorem, 
FT. Surprisingly perhaps, it is valid far from equilibrium in the nonlinear response regime 
[1]. Since 1993 there have been a number of derivations and generalisations of the FT. 
Evans and Searles [3-5] gave a derivation, similar to that given here, which considered 
transient, rather than steady state, nonequilibrium averages and employed the Liouville 
measure. Gallavotti and Cohen [6,7], gave a proof of the formula for a nonequilibrium 
stationary state, based on a Chaotic Hypothesis and employing the SRB measure. In the long 
time limit, when steady state averages are independent of the initial phase used to generate the 
steady state trajectory, averages over transient segments which originate from the initial 
equilibrium microcanonical ensemble can be expected to approach those taken over 
nonequilibrium steady state segments. Thus for chaotic systems both approaches should be 
able to explain the steady state results. However this point is being debated [8]. Other 
generalisations of the FT have recently been developed [9-12]. 

In a footnote to their original paper ECM2 also pointed out that in the weak field regime, 
there was a connection between the FT, the Central Limit Theorem (CLT) [13, 14], and 
Green-Kubo relations [1, 4]. In the present paper we explore this connection further and 
consider the validity of the Green-Kubo relations far from equilibrium. We show that if the 
distribution of the time averaged dissipative flux, p( J(t)) , is Gaussian arbitrarily far from the 
mean, then from the FT one can derive both generalised Einstein and generalised Green-Kubo 
relations for the relevant transport coefficient. Both isothermal (ie isokinetic) and 
isoenergetic dynamics are considered. We conduct computer simulations which prove that 
outside the linear regime, these generalised Green-Kubo and Einstein relations are incorrect. 



It turns out, that in order for the nonlinear Green-Kubo relations to be valid p(J(t)/Oj ) must 
be a normalised Gaussian when both t and J(t)/aj -^ oo . However, this is not guaranteed 
by the Central Limit theorem [15] and the nonlinear Green-Kubo relations are invalid. 



II. NEMD DYNAMICAL SYSTEMS 

The development of NonEquilibrium Molecular Dynamics, NEMD, over the previous two 
decades has lead to a set of deterministic algorithms (ie N-particle dynamical systems) from 
which one can in principle, calculate correct values for each of the Navier-Stokes transport 
coefficients [16]. These dynamical systems actually duplicate the salient features of real 
experimental nonequilibrium steady states. In the linear regime close to equilibrium, 
nonequilibrium statistical mechanics is used to prove that in the large system limit the 
calculated transport properties are correct. Using NEMD one can calculate far more than just 
transport coefficients. One can also correctly calculate the changes to the local molecular 
structure and dynamics, caused by the applied external fields. 

Consider an N-particle system in 3 Cartesian dimensions, with coordinates and peculiar 
momenta, {qj,q2,..qN,Pi,..PN} = (q,p) = F. The internal energy of the system is 

N 

Hq = V Pi / 2m -l-$(q) where 4>(q) is the interparticle potential energy which is a function 

i=l 

of the coordinates of all of the particles, q. In the presence of an external field Eg, the 
thermostatted equations of motion are taken to be, 

qi=Pi/m + Ci(r)E, 

Pi = Fi(q) + Di(r)Ee-a(r)Pi (1) 

where 

Fi(q) = -a4>(q)/aqi (2) 

and a is the thermostat multiplier derived from Gauss' Principle of Least Constraint in order 

N 

to fix the peculiar kinetic energy, K = V Pj / 2m, or the internal energy, Hq. In a constant 

i=l 

energy system the thermostat multiplier is easily seen to be, 

aE=-J(r)VE,/2K (3) 

while in an isokinetic system the corresponding expression for the multiplier is, 



£P^.(Fi+D,FJ 
«K= "^ /2K- (4) 

We note that the thermostatted equations of motion are time reversible. The dissipative flux 
is defined in terms of the adiabatic (ie unthermostatted) derivative of the internal energy, 

<^-J(r)VFe (5) 

where V is the system volume. 

In an isokinetic system, the balance between the work done on the system by the external 

field and the heat removed by the thermostat implies that 

t t 

lim fdsJ(r(s))VFg=-2Ko lim fdsaK(r(s)), (6) 

while in a constant energy system, energy balance is exact instantaneously, 

J(r)VF, = -2K(p)aE(r). (7) 

In equation (6), Kq is the (fixed) peculiar kinetic energy, 

3NkBT/2 = 3N(3oV2 = Ko- (8) 

A shorthand notation will be used to refer to the time averaged value of a phase function 

along a trajectory segment, r+(s); < s < t. We will write, 

t 
A+(t)^)/jdsA(r+(s)). (9) 



Since the dynamics is time reversible, for every trajectory segment r_,_(s); < s < t, there 
exists an antisegment, r_(s);0<s<t. A plus or minus sign is ascribed to a particular 
trajectory segment depending on the sign of the time averaged value of the thermostat 
multiplier: therefore by definition a_,_(t) > 0. The time reversed conjugate of the segment 
r_,_(s); < s < t, namely r_(s); < s < t, is termed an antisegment and. 



A_(t)^)/jdsA(r_(s)). 



(10) 



Depending on the parity of the phase function A(r) under the time reversal mapping there 
may be a simple relation between A_,_(t) and A_(t). Without loss of generality we take the 
external field to be even under time reversal symmetry, therefore the dissipative flux is odd 
and. 



J_(t) = -J+(t),Vt. 



(11) 



Using this notation the dissipative flux is related to the phase space compression 
accomplished by the thermostat. 



lim |3J+(t)VFg = - lim 3Na+(t) 

(t^oo) (t->oo) 



pJ (t)VFe=-3Na+(t) 



isokinetic 



isoenergetic 



(12) 



where for both the isokinetic and isoenergetic systems, 
|3J(r)V = 3NJ(r)V/2K(p) 



but in the isokinetic case, the peculiar kinetic energy K is a constant of the motion. Since (3 is 
always positive, we see from (12), that the sign convention for distinguishing segments and 
antisegments can equally well be taken from the sign of the dissipative flux. 



III. THE (TRANSIENT) FLUCTUATION THEOREM 

For our system, since the adiabatic incompressibility of phase space (AIF) tBlds [16], the 
Liouville equation for the N-particle distribution function f(r,t) , reads, 

^^^ = -f(r,t)5r •/ = 3Na(r)f(r,t) + o(i)f(r,t). (i3) 

dt / "^ 

The 0(1) terms are omitted in the following discussion. Incorporation of these terms poses 
no difficulty but complicates the expressions and the consequences can be neglected in the 
large system limit. The solution of this equation can be written as [4] 

t 
f(r+(t),t) = exp[J3Na(r+(s))ds]f(r+,0) 



(14) 
= exp[3Na+(t)t]f(r+,0) 

This is known as the Lagrangian form of the Kawasaki distribution function [4]. 

Consider the propagation of a phase point along a trajectory in phase space. If we select 
an initial, t = 0, phase, r(i), and we advance time from to x using the equations of motion (1) 
we obtain r(2) = r(x;r(i)) = exp[iL(r(i),Fg)x]r(i), where the phase Liouvillean, iL(r(i),Fg), is 
defined as, iL(r,FJ = q(r,FJ •3/3q + p(r,FJ«3/3p. Continuing on to 2x gives r(3) = 
exp[IL(r(2),FJJx]r(2)= exp[IL(r(i),Fg)2x]r(i). This is demonstrated in Figure 1. 

From this trajectory segment, we can construct a time-reversed trajectory. At the midpoint 
of the trajectory segment r(i 3) (i.e. at t = x) we apply the time reversal mapping, M(T), tB" r(2) 
generating M(T)r(2)& ^"(5^ If we now propagate backward in time keeping the same external 
field, we obtain r(4) = exp[-irXr(5),Fg)x]r(5). r(4) is the initialIT'= phase from which a 
segment r(4 6) can be generated with r(5) = exp[iE,(r(4),Fg)2x]r(4). We denote the trajectory x- 
segment r(i)^r(3^ as r^3)=r_,_. similarly r^4g^=r_. Using the symmetry of the 
equations of motion it is trivial to show that, I(r(2)) = -J(r(5)) and that J(t; r_,_, < t < 2x) = 
-J(2x-t; r_, < t <2x), - see Figure 1. We now have an algorithm for finding initial phases 
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which will subsequently generate the conjugate segments. 

The ratio of probabilities of finding the initial phases r(i), r(4) which gEhSfate these 
conjugate segments will now be discussed. In a causal universe, the probabilities of 
observing the segments r+ and P. are proportional to the probabilities of observing the initial 
phases which generate those segments [3, 5]. It is convenient to consider a small phase space 
volume, 5V(r(i)(0)) about an initial phase, r(i)(0). If we are conlidering isoenergetic 
dynamics, then the initial equilibrium phases are distributed microcanonically, and therefore 
the probability of observing ensemble members inside 5V(r(i)(0)), is proportional toF 
5V(r(i)(0)) (ffir generalisations to other ensembles see [17]). From the Liouville equation 
(13) and the fact that for sufficiently small volumes, 5V(r(t))F l/f(r(t),t),Iwe can make the 
following observations: 5V2 = 5Vi(x) = 5Vi(0)exp[-f^3Na(s;rn))cls] and, 5V3=5Vi(2x)= 
oVi(0)exp[- 3Na(s;r/iOds] . Because the segment r(^6) is related to r(il3) by MT which is 
applied at t = x, and the Jacobian of MT is unity, 5V2 = SV5 ^> 5V3 = 5V4 and 5Vi(0) = 5V6. 

However, since 5Vi(0) and 5V4(0) are volumes at t = and since the distribution of initial 
phases is microcanonical, we can compute the ratio of probabilities of observing t = phases 
within 5Vi(0) and 5V4(0). This ratio is just the volume ratio, 

5V4(0)/5Vi(0) = 5Vi(2x)/5Vi(0) = exp[j^-3Na(s;r(i))ds], Vx. (15) 

There may be trajectory segments whose initial phases lie outside the phase space volume 
5Vi(0) but which have the same value of a(2t) as those lying inside 5Vi(0). Suppose that 
there are two noncontiguous subvolumes of phase space 5Vi(0), 5Vi'(0) from which 
trajectories originate which, after a time 2x , have time averaged values of a which lie in the 
range: (a^_(2x),a^_(2x) + da), (see Figure 2). Suppose that at time 2x these volumes evolve 
to: 5Vi(2x)=5V3, 5Vi'(2x)=5V3'. At time zero the corresponding volumes for the 
corresponding antitrajectories are: 5Vi(2x)=5V4, 5Vi'(2x)=5V4'. 

The ratio of probabilities of observing trajectories (a^.(2x),a_,_(2x) + da) compared to the 
corresponding antitrajectories, (a_(2x),a_(2x) + da) , is 



(5V4 + 6V4. )/ _ (6V3 + 6V3. )/ 

/(6Vi + 6V1. ) ~ /(6V1 + 6V1. ) 

(6V1 exp[3Na+(2x)2x] + 6V1. exp[3Na+(2x)2x]) / 
lim da^o /(SVi + 6V1. ) 

= exp[3Na+(2x)2x] 

(16) 

since a(2x) is the same for both trajectories. This shows that even when noncontiguous 

regions of phase space have the same time averaged values for the thermostatting multiplier, 

the ratio of probability that a(t) = a_,_(t) = A to the probability that a(t) = a_(t) = -A, (ie 

p(a(t) = A)/p(a(t) = -A), where A is any required value of the time average of a is, 

P(^('^ = ^^ =exp[3NAt] (17) 

p(a(t) = -A) ^^ J- ^ '^ 

This formula is exact for transient trajectory segments of a system undergoing isoenergetic 

dynamics. For isoenergetic dynamics there is a linear relationship between the value of a and 

the value of (3J, so 

p(a(.) = A) ^ p(jI(.) = B) , pf,A,] = exp[-BVF,.]. (18) 

P(a(t) = -A) p(pj(t) = -B) " "^^ " 

where B=-3NA/VFe. For isokinetic dynamics, the procedure above can be used to show that 
[17], 



Pa^^^^ = exp[-BpVF,t]. (19) 

p(J(t) = B) ^ ^^ ^ ' 



If we are interested in steady state segments, equations (18) and (19) will only be true in the 
limit as the segment duration, t, goes to infinity [1, 6-8, 17] and only when the steady state is 
unique: 



lim -In 

(t^c«)t 



^p(|3J(t) = B) 



VJ 



= -BVF,. (20) 



p(|3J(t) = B)^ 
Equations (18, 19 and 20) express what has become known as the Fluctuation Theorem, (FT), 
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for the dissipative flux for isokinetic and the isoenergetic dynamics [1, 2-7, 17]. 



11 



IV. EINSTEIN AND GREEN-KUBO RELATIONS 

We consider first the isokinetic case. In this case P is a constant of the motion: 
(3J(t) = PoJCt). It might be expected that as the averaging time, t becomes arbitrarily large 
compared to the Maxwell time, Xj^, which characterises serial correlations in the dissipative 
flux, contributions to the trajectory segment averages of the dissipative flux, { J(t)}, would 
become statistically independent and therefore satisfy the Central Limit Theorem, (CLT). 
That is, as t^oo, the distribution would approach a Gaussian. If the distribution is 
Gaussian, it is trivial to show that there is a relation between the logarithm of conjugate 
probabilities of time averaged steady state dissipative fluxes and the variance of the 
distribution of those averaged dissipative fluxes, 



lim -In 

(t^oo)t 



r „/i^/.^ OT^^ A ^ A„/T/.^ Tl^^ 



p(PJ(t) = PB) 
p(|3J(t) = -pB) 



= lim -In 

(t^oo) t 



p(J(t) = B) 
p(J(t) = B) 



(21) 

2B(J)p 

= lim — ^ 

(t--)taj(^/t) 

2 — 

where Ot, ,(t) is the variance of the distribution of { J(t)}. Combining this equation with 

J\\.) 

(20) shows that if the distribution is Gaussian there must be a trivial relation between the 
variance and the mean of the distribution of averaged fluxes [18]. From this relation the 
nonlinear transport coefficient is given , 

L(Fe) = -—^ = lim iPoVta? (22) 

Fg (t^oo)^ •'*-'••' 

In the zero field limit this equation constitutes an Einstein relation for the linear transport 
coefficient, L(0). Except for the case of colour conductivity where (22) is equivalent to the 
standard Einstein expression for the self diffusion coefficient [19], these zero field Einstein 
relations are not well known. For nonzero applied fields, the generalised Einstein relation for 
the field dependent transport coefficient, h(F^), (22) is, as we shall see, incorrect. 

In the long time limit the variance of the steady state distribution of t-averaged fluxes. 
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4t)(Fe) = ((J(t)-(J>F„)')p , (23) 

satisfies the equation (see [4] and also the Appendix), 

lim<,(F.) = 2L,(0;F.)/ .^2Li(0;F.)/ 

(t^oo) ■>('-> /PO^ (t^oo) /PO 



.vt 



(24) 



_ 2Lj(0;FJ 



where 



Lj(s;F,) ^ poV jdt e-^'((J(0) - (J)p^ )(J(t) - (J)p^ )). (25) 



Lj(s;Fg) is the frequency and field dependent Green-Kubo transform (GK), of the dissipative 
flux and therefore is essentially the Fourier-Laplace transform of the field dependent 
dissipative flux autocorrelation function evaluated in a NESS with an applied field Fg. 

L;(s;FJ^^Lj(''Pe^3 . (26) 

The factor po is included in the GK transform to make the expression consistent with the 
Green-Kubo expression for the transport coefficient at zero applied field. 

Combining (22) and (24), shows that if the t-averaged dissipative fluxes are Gaussian, then 
the non/mear phenomenological transport coefficient, L(Fg), is given by the zero frequency 
Green-Kubo transform of the dissipative flux, 

oo 

L(FJ = Lj(0;FJ = PqV jdt ((J(0) - (J)p^ )(J(t) - (J)p^ ))^ . (27) 

In the zero field limit (27) reduces to the correct well known Green-Kubo expression for the 
linear transport coefficient, L(0). The relationship between the FT and GK expressions in the 
linear regime has been considered previously [4, 10, 20-22]. In the present paper, simulations 
are carried out to test these relationships in the nonlinear, large field regime. These numerical 
calculations show that this generalised Green-Kubo relation (27), is not valid, forcing us to 
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conclude that the distribution is not sufficiently Gaussian far from equilibrium and far from 
the mean. 

In the isoenergetic case, if the distribution is Gaussian, we have. 



lim -In 

(t^oo)t 



/ 



P(jj(t) = B) 
p(|3J_(t) = B) 



2B(|3J). 



lim 



pj(t) 



lim 

(t^oo) to 



-2B(|3J)j 



(28) 



pJ(t) 



Combining this equation with (20) shows that if the distribution is Gaussian there there must 
again be a trivial relation between the variance of the distribution of averaged fluxes and the 
nonlinear transport coefficient. 



-(PJ)f 

(P)p L(FJ ^ — — ^ = lim iVto^ 



(t^cx.) • 



pj(t) 



(29) 



Were such a relation to be true at large fields it would constitute a generalised Einstein 
relation for the field dependent transport coefficient L(FJ . In the long time limit the variance 
of the steady state distribution of t-averaged fluxes. 



cyL/FJ = ((|3J(t)-(|3J)p)'^^ 



pJ(t)' 



(30) 



satisfies the equation 



(t^oo) 



.2_ 
*pj(t)^^e> 



2(|3)p Lj(0;FJ 



lim to^^ rFJ= ' '^^ - Ar + lim 



2(P)p L'j(0;FJ 



(t^oo) 



2(|3L Lj(0;Fe) 



where 



(31) 



(|3)p Lj(s;FJ ^ VJdt e-^V(|3J(0)-(|3J)p^)(|3J(t)-(|3J)p ) 



(32) 



Combining (28) and (30), shows that if the distribution of the t-averaged dissipative fluxes is 
Gaussian, then the nonlinear phenomenological transport coefficient, L(Fg), is given by the 
zero frequency Green-Kubo transform of the dissipative flux, 
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L(F,) = Lj(0;FJ ^ V(|3);^' jdt ((|3J(0)-(|3J)p^)(|3J(t)-(|3J)p^)). (33) 

' 

Not surprisingly, results of numerical tests of this relationship indicate that (33) is also not 
correct. 
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V. NUMERICAL RESULTS 

Steady state NEMD simulations of a fluid undergoing shear flow were used to test the 
accuracy of the expressions derived above. All simulations were carried out in two Cartesian 
dimensions with interactions between particles given by the Weeks-Chandler-Anderson 
repulsive pair potential. Note that Lennard-Jones reduced units are used in the figures and 
throughout this section. In both cases, simulations were carried out for systems of 200 
particles and for the isokinetic system, the temperature was constrained at T = 1 .0, whereas 
for the isoenergetic system the internal energy was constrained at E/N = 1.56032. For the 
isokinetic fluid, two densities, n=N/V, were considered: n = 0.4 and n = 0.8; and for the 
isoenergetic fluid the density was set to n = 0.8. 

The SLLOD equations of motion with Lees-Edwards periodic boundary conditions were 
employed to model the shear flow, and a Gaussian thermostat or ergostat used to maintain a 
steady state [16]. The adiabatic SLLOD equations give an exact representation of shear flow 
arbitrarily far from equilibrium and Lees-Edwards periodic boundary conditions give the 
unique generalisation of periodic boundary conditions to planar Couette flow. The SLLOD 
equations (analogous to equation (1)) are given by: 

Pi=Fi-iYPyi-api 

where y is the strain rate and a is the isokinetic or isoenergetic thermostat multiplier. When 
the kinetic energy is a constant of motion, 

N 

E^i-P.-TP^Pyi 

Wk = "^^-^ (35) 

£pi Pi 

i=l 

while if the internal energy is a constant of motion. 
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-yP V 

«E = N "' (36) 

Lp. Pi 

where Pxy is the xy element of the pressure tensor, 

N N 

Pxy V = £ PxiPyi - i £ XijFyij , (37) 

i=l i,j=l 

which is the dissipative flux: J = Pxy. The nonlinear shear viscosity, ri(Y) is the nonlinear 
transport coefficient calculated using this algorithm. We note that in contrast to the 
discussion above, the dissipative flux for shear flow is even under the time reversal mapping, 
MT(x, y, Px, Py) = (x, y, - Px, - Py) and the strain rate is odd. However we can choose the 
strain rate to be even and the dissipative flux odd by employing the Kawasaki mapping [16], 
MK(x, y, Px, Py) = (x, -y, - Px, Py). 

Firstly we carried out simulations to show that in the long time limit, for an isokinetic 
system the variance of the distribution of { J(t) } is related to the zero frequency Green-Kubo 
transform of J(t) by equation (24), and for an isoenergetic system, the variance of the 
distribution of { J(3(t) } is related to the zero frequency Green-Kubo transform of J|3(t) by 
equation (31). The behaviour at various strain rates was examined and the results are shown 
in Figure 3. Equations (24) and (31) are found to be verified in all cases. 

We tested the nonlinear Green-Kubo relations (27) and (33) for the systems described 
above and the results are shown in Figure 4. Since L (Fg) is also related to the variance of the 
distribution of {J(t)} or {J|3(t)} as t ^ oo for the isokinetic or isoenergetic system, 
respectively, results obtained directly from the variance are also presented. Clearly the 
equivalence of L(Fe) and L(Fe) is only observed at small fields. At intermediate fields the 
Green-Kubo transform of the dissipative flux L(Fe), underestimates the actual transport 
coefficient while at high fields L(Fg), overestimates the transport coefficient. We conclude 
that nonlinear Green-Kubo relations (27, 33) are not valid in the far from equilibrium regime. 
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VI. DISCUSSION 

In the zero field limit, thermostatted linear response theory can be used to determine the 
field dependent transport coefficients. For the isokinetic response: 

lim L(FJ = - lim ^ = |3„V jdt(J(0)J(t))„^ (38) 

where the ensemble average { )q^ is over the equilibrium isokinetic ensemble. This 
expression derived from linear response theory is identical to (27) in the limit F^^O, which 
was derived using both the CLT and the FT. The results in Figure 4 confirm the agreement of 
(38) and (27) in the zero field limit with linear response theory. 

This work also shows that in the zero field limit, one can calculate linear transport 
coefficients by considering the limiting long time variance, Oj (F^ =0), of the distributions 
of J(t) (24, 31), rather than by computing autocorrelation functions of the dissipative flux and 
then performing the appropriate long time integrals. The variance of the t-averaged flux 
therefore provides an alternative route to the linear transport coefficients and equations (24, 
31) thus provide useful Einstein routes to linear transport coefficients. 

We now turn to the question of why the nonlinear Green-Kubo and Einstein expressions 
fail, far from equilibrium. A necessary condition for the CLT is the statistical independence 
of the sample averages. A breakdown of this independence could be responsible for the 
breakdown of the CLT, Einstein and Green-Kubo expressions. However, trajectory segments 
that are much longer than the Maxwell time, Xm, which characterises the decay of the 
autocorrelation function of the dissipative flux autocorrelation function, should have no 
correlations between successive samples of J(t). Thus, regardless of the distribution of J(t) 
we expect that for long enough t, the CLT will apply. Figure 5 compares the decay time 
autocorrelation function of J(t) for different applied fields. At moderate fields Xm is less than 
it is at equilibrium. Only at very large fields does Xm increase. This means that possible 
decay time divergences or anomalies are not responsible for the breakdown of the nonlinear 
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Einstein and GK expressions (33, 27). Further, if one computes the distribution of J(t), for 
various values of t, one cannot observe departures from Gaussian behaviour for values of t » 
Xm, in the neighbourhood of the mean current. 

Figure 6 compares the distributions of J(t) (the distribution of the instantaneous flux) and 
J(t) for an equilibrium and nonequilibrium system with a strain rate which ensures it is in the 
nonlinear regime (T=1.0, n=0.8, '^1.0). While the skewness, yi, and kurtosis, K, for the 
instantaneous equilibrium J(t) distribution are zero within error bars which is consistent with a 
Gaussian distribution, the skewness is non-zero for the sheared system (yi = -0.23+0.01, K = 
0.13+0.04, respectively). The distributions of the time-averaged fluxes, J(t), were obtained 
for a trajectory segment of length t = 4.0 and both distributions, as expected, appear 
Gaussian. The skewness of the distribution for the sheared system is yi = -0.064+0.004, and 
the kurtosis K = -0.02+0.02. Thus on the basis of these tests although for a sheared system 
the distribution of J(t) is not Gaussian the distribution of J(t), for a trajectory segment of 
length t = 4.0, is on the scales shown in Figure 6, already indistinguishable from a Gaussian. 

As noted in references [21, 22], the distribution of J(t) and J(3(t) cannot be exactly 
Gaussian because the values of these variables are bounded. In practice however these 
bounds are so large that they become irrelevant in the limit t ^ oo where the t- averaged 
distributions collapse to zero variance distributions. Moreover, the bounds still apply in the 
zero field limit where the Green-Kubo and Einstein expressions are all valid. Thus the 
boundedness of the fluxes cannot be the responsible for the breakdown of the nonlinear GK 
and Einstein expressions. 

If we examine the derivation of equations (27) and (33) more closely, it can be seen that in 
order to obtain a GK expression we require the distribution at both J(t) = J_,_(t) and 
J(t) = -J_(t) be well approximated by a Gaussian for times sufficiently long that the GK 
integrals have converged, t » XM(Fe) (see Appendix for details) [13, 14]. Any deviations 
from the behaviour indicated by (21) and (28) will be related to the relative deviation of the 
distribution from a Gaussian at both J_,_(t) and J_(t) . It is therefore of interest to consider the 
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rate of convergence to a Gaussian. The magnitude of the relative deviation of the distribution 
p((J(t) - JyoT, ) from a normalised Gaussian generally increases with the separation of 
J(t) from the mean J for sufficiently large separations (see, for example, section 7.2 of [13]). 
Here Oy, is the standard deviation of the distribution p(J(t)) when t = Xm- In the t ^ oo 
limit, at fixed J(t), the magnitude of the relative devaition of p(J(t)) from a Gaussian 
becomes infinite. This means that in the t ^ oo limit, the CLT gives information which is not 
sufficiently precise to derive Green-Kubo relations for non-zero applied fields. 

We illustrate this point in more detail. Suppose that J(t) = J+(t) is equal to the mean 
current, J; then the conjugate trajectories will have J(t) = J_(t)= -J. Clearly 
I J_(t) - J 1= 2L(FJF^ . Therefore using equation (45) of Appendix A, we find in the t ^ oo 
limit, except when F^ = , 



J_(t)-J 



/ Oj, , = 2L(FJFe / Oj, , - FeV2VtL(FJ|3 ^ - . (39) 



For any non zero field, if J(t)= J, then as t increases, the value of J_(t) moves further and 
further into the wings of the normalised distribution where the magnitude of the relative 
deviation of p(J(t)) from a Gaussian grows without bound. Strictly speaking therefore, in the 
infinite time limit, for any finite field, the relative deviation of p(J(t)) from a Gaussian, 
evaluated in the neighbourhood of the mean anticurrent, -J grows without bound and 
nonlinear Green-Kubo relations cannot be derived. However, in practice one does not need to 
take the infinite time limit. Considering the shift in the mean value of the dissipative flux 
with field shows that the nonlinear GK expression will be approximately correct provided, 

Fe < o(l / V|3VMXM(Fe)L(Fj), (40) 

where V^ is the minimum volume required for transport coefficient to be approximately 
equal to its large system, limiting value [24]. Clearly the nonlinear GK relations satisfy this 
relation only in a small neighbourhood including F^. = . For the systems studied here, 
equation (40) predicts that the nonlinear GK relations will be approximately correct provided 
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Y <~ 10 '. This is in agreement witli experimental observations given in Figures 4(a), (b), 
4(c). 
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Appendix. 

The variance of the time- averaged dissipative flux is given by, 



4t) = ((J(t)-(J)r 

= ^^(£'ds,AJ(s,))(£ds,AJ(s,))^ (40) 



where AJ(t) = J(t) - (J) . Using a change of variables: Xj = Sj - s^ and Xj = Sj + Sj this 
integral can be written: 

4t) = Wat^^^r ds,(AJ(i(x, + x,))AJ(i(x, - X,))) 

(41) 



+ — j^dx^^ ds2(AJ(i(Xi + X2 - X))AJ(|(X2 - X; - x))) 



2t^Jo 

Since correlation functions are invariant under a time translation in the steady state, and 
using the symmetry of the functions we obtain, 

''lr^fo^^2iy{AKx,)Am). (42) 

Changing the order of integration gives: 

2 r' 



4o = elo^'J/'^^^J^''^^J^^^) 

= -^j^'dXi^ dX2(AJ(x,)AJ(0)) . (43) 

= ^rdXi(AJ(x,)AJ(0)(t-Xi)) 
t •^'' 



Therefore, for any steady state system, at all times: 

to^^^,(FJ = 2 jds((J(0) - (J)p^ )(J(s) - (J)p^ )} - ^ jds((J(0) - (J)p^ )(J(s) - (J)p^ )}s. (44) 

At any time greater than the time required for the time correlation function to decay to zero, 
tc > tM, in an isokinetic system, f ^ ds/(J(0) - (J)p )(J(s) - (J)p )\ = Lj(0;FJ / (P^V) and 
fj ds((J(0) - (J)p^ )(J(s) - (J)p^ )}s = -L;(s;FJ / (PoV). Therefore, 
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tc<)(I^)^ 



2Lj(0;Fj/ ,2L;(0;FJ 



PoVtc • 



(45) 



7/' the distribution is Gaussian at J(t) and - J(t) at tc, then assuming that the second term of 
(45) is negligible and that the FT is true at t = tc, combining (20) and (45) gives, 



-(j: 



~ 2 



PoVtcO;,^,=L,(0;FJ. 



(46) 



That is, a GK expression is valid. 
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Figure Captions 

Figure 1. The shear stress, Pxy for trajectory segments from a simulation of 200 disks at T = 
T. Zpi2/2mNkB = 1.0, and n = N/V = 0.4. The trajectory segment, r(i 3), wis obtained from a 
forward time simulation. At t = 2, a time reversal map was applied to r(2), tFgive r(5) (Tor 
the SLLOD equations of motion the time reversal map is the Kawasaki map (x, y, px, Pz, Y)^ 
(x, -y, -px, Pz, Y)). Forward and reverse time simulations from this point give the trajectory 
segmentJT(5 6) and r(5 4), respectively. If one inverts Pxy in Pxy = and inverts time about t = 
2, one transforms the Pxy(t) values for the anti-segment r(4 5) into tlose for the conjugate 
segmeiT, r(i3). 

Figure 2. Schematic diagram showing two disconnected subvolumes of phase space 5Vi(0), 
5Vi'(0) from which trajectories originate which, after a time 2x , have time averaged values of 
a which lie in the range: (a_,_(2x),a_,_(2x) + da) . At time 2x these volumes evolve to: 5Vi(2x) 
= 5V3, 5Vi'(2x) = 5V3'. See equation (16). 

Figure 3. Calculation of Lj(0;Fe) from the variance of distributions of the t-averaged 
dissipative flux for simulations at: (a) constant temperature with T = 1.0, n = 0.4, Fe = 0.0 
(unfilled circles), Fg = 1.0 (filled circles) and Fe = 2.0 (squares); (b) constant temperature 
with T = 1.0, n = 0.8; Fe = 0.0 (unfilled circles), Fe = 0.5 (filled cirlces) and Fe = 1.0 
(squares); and (c) constant internal energy with E/N = 1.56032, n = 0.8, Fg = 0.0 (unfilled 
circles), Fg = 0.3 (triangles) and 0.5 (filled circles). The crosses show Lj(0;Fe) determined 
from the zero frequency Green-Kubo transform of the dissipative flux. 

Figure 4. The filled circles show the viscosity as a function of strain rate for systems at (a) 
constant temperature with T = 1.0 and n = 0.4; (b) constant temperature with T = 1.0 and 
n = 0.8; and (c) constant internal energy with E/N = 1.56032 and n = 0.8. The crosses are 
predictions determined from the Green-Kubo expression (equation (27) for the isokinetic case 
and equation (33) for the isoenergetic case) and the squares are predictions from the variance 
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(equation (22) for the isokinetic case and equation (29) for the isoenergetic case). 

t 
Figure 5. The decay of the shear stress time correlation function, J ds^P^ (0)Pxy(s)V at 



equilibrium and in nonequilibrium steady states for a systems at constant temperature with T 
= 1.0 and n = 0.4. The full line is at equilibrium (y= 0) the dotted line is for a simulation with 
a strain rate of y = 1 and dashed line with a strain rate of y = 2. 

Figure 6. (a) The instantaneous (small circles) and time-averaged distribution (squares, t = 
4.0) of the dissipative flux for an equilibrium system at T = 1.0 and n = 0.8. (b) The 
instantaneous (small circles) and time-averaged distribution (squares, t = 4.0) of the 
dissipative flux for a nonequilibrium system at T = 1.0, n = 0.8 and with y= 1.0. 
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